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This  technical  report  describes  computer  code  generated  for  the 
study  of  dipole  and  finite  element  models  of  magnetic  field/defect 
interactions  associated  with  leakage  field  methods  of  nondestructive 
testing.  It  is  the  principal  Investigator's  contention  that  finite 
element  analysis  techniques  could  be  used  to  develop  defect  characteri- 
zation schemes  for  all  electromagnetic  methods  of  nondestructive  testing, 
as  well  as  providing  valuable  insight  into  the  study  of  defect  signal 
phenomena . 

Experimental  results  are  given,  showing  Hall  probe  generated  leak- 
age field  profiles  around  rectangular  slots  in  a steel  bar,  and  compared 
with  both  dipole  and  finite  element  model  predictions.  The  finite  ele- 
ment code  has  already  been  used  to  produce  a library  of  leakage  field 
profiles  for  a variety  of  surface  and  subsurface  defects  in  a ferromag- 
netic bar  c::rrying  an  axial  dc  magnetizing  current  and  as  part  of  ARO 
contract  DAAG29-76-G-0249,  the  code  will  be  extended  to  Include  residual 
forms  of  magnetization. 


L 


1 


Introduction 

During  the  past  decade,  the  urgent  need  for  Incorporating  nondestruc- 
tive testing  and  evaluation  considerations  directly  Into  the  engineering 
design  process  has  been  amply  Illustrated  by  material  failures  In  aero- 
space, electric  pcwer  eind  transportation  Industry  equipment.  The  economic 
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Impact  of  such  failures  Is  well  documented  and  provides  a major  Impetus 
to  Improve  all  aspects  of  the  nondestructive  testing  art.  Considerable 
progress  has  been  made  toward  this  end  through  such  efforts  as  the  ASPA/ 
AFML  program^.  Although  the  work  has  concentrated  primarily  on  ultrasonic 
techniques,  much  of  the  research  philosophy  developed  for  the  program 
with  regard  to  the  study  of  basic  phenomena^,  development  of  models^, 
signature  Identification  by  signal  processing  and  the  subsequent  accept/ 
reject  decision  founded  on  a knowledge  of  fracture  mechanics  and  related 
failure  probabilities  could  and  should  be  applied  to  other  nondestructive 
testing  techniques.  The  cornerstone  of  such  an  approach  Is  the  develop- 
ment of  an  adequate  mathematical  model  for  the  study  of  the  basic  field/ 
defect  Interactions.  Such  a model  Is  needed  In  order  to  develop  a defect 
characterization  scheme  and  to  Identify  suitable  parameters  for  the  signal 
processing. 

Despite  recent  developments  In  automatic  defect  characterization 

associated  with  eddy  current  and  leakage  flux  methods  of  nondestructive 
9 10 

testing  ’ , the  subject  of  electromagnetic  methods  of  nondestructlvely 

testing  ferromagnetic  materials  Is  characterized  largely  by  empirical 
knowledge.  Where  closed  form  mathematical  solutions  do  exist,  describing 
electromagnetic  fleld/defect  Interactions,  the  underlying  assumptions  of 
the  theories  tend  to  Invalidate  any  realistic  application  of  the  results 
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to  the  problem  of  defect  characterization  in  ferromagnetic  materials. 

The  problem  of  modeling  electromagnetic  fleld/defect  interactions 
in  ferromagnetic  materials  is  complicated  not  only  by  the  nonlinear 
magnetization  characteristic  of  the  material  and  the  awkward  defect 
boundaries  but  also  by  the  number  of  different  phenomena  utilized^^. 

Figure  1 shows  the  three  major  modes  of  field/defect  interaction.  In 
magnetic  particle  and  magneto graphic  methods  of  nondestructive  testing, 
the  ferromagnetic  specimen  to  be  tested  is  initially  magnetized  with  dc 
current.  Residual  leakage  fields  are  set  up  around  defects  present  in 
the  specimen.  The  presence  of  such  fields  can  be  observed  either  direct- 
ly using  magnetic  particles  or  indirectly  by  recording  the  leakage  fields 
on  magnetic  tape.  With  reluctance  and  perturbation  techniques,  a constant 
dc  magnetization  current  is  used  to  set  up  active  leakage  fields  around 
defects  in  the  parts  undergoing  Inspection,  which  can  be  detected  using 
any  flux  sensitive  transducer^^.  Finally  in  eddy  current  testing,  ac 
excitation  is  used  to  induce  secondary  currents  and  fields  in  the  specimen. 

In  this  case  the  presence  of  a defect  causes  changes  in  both  induced  cur- 
rents and  fields,  resulting  in  measurable  Impedance  variations  in  the  pick- 
up coll. 

Differences  in  the  residual,  active  and  eddy  phenomena  are  summarized 
in  Figure  1 which  shows  the  regions  of  the  material's  B-H  characteristic 
over  which  the  various  phenomena  occur. 

The  major  purpose  of  this  technical  report  is  to  give  details  of  a 
finite  element  model  for  predicting  leakage  fields  around  defects  in 
ferromagnetic  parts  carrying  dc  magnetization  current.  It  is  the  prin- 
cipal Investigators  contention  that  the  finite  element  model  could  be 
used  as  the  basis  for  developing  defect  characterization  schemes  for  all 

I 

...  J 


3 

electromagnetic  methods  of  nondestructive  testing.  Including  active 
direct  current,  residual  and  alternating  current  methods  of  excitation. 
The  following  sections  describe: 

1.  Actual  residual  and  active  leakage  field  profiles  for 
rectangular  slots  in  a steel  bar  obtained  by  scanning 
a Hall  probe  over  the  surface  of  the  bar. 

2.  Dipole  model  predictions  for  active  and  residual  leakage 
fields  around  the  same  slots. 

3.  Finite  element  predictions  for  the  active  leakage  fields 
around  the  same  slots. 

report  concludes  with  a discussion  of  the  relative  merits  of 
the  two  modeling  techniques  and  indicates  how  the  finite  element  model 
might  be  extended  to  include  residual  magnetization  and  eddy  current 
phenomena.  An  extensive  bibliography  is  included  for  further  reference. 

Active  and  Residual  Leakage  Fields  Around  Rectangular  Slots 

The  general  shapes  of  residual  and  active  leakage  fields  around 

defects  in  ferromagnetic  materials  are  well  known  and  have  been  widely 

12“16 

reported  in  the  literature  . However,  in  order  to  provide  a basis 
of  comparison  with  theoretical  modeling  studies,  residual  and  active 
leakage  field  profiles  were  monitored  aroimd  rectangular  slots  in  a 
steel  bar  using  the  apparatus  shown  in  Figure  2.  A miniature  Hall  probe 
was  used  to  measure  the  vertical  component  of  the  leakage  field  flux 
density.  Leakage  field  profiles  were  obtained  by  driving  the  probe  with 
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a uniform  velocity  over  the  surface  of  the  bar  using  an  x-y  positioning 
mechanism.  The  output  of  the  probe  was  fed  to  a Gaussmeter  having  an 
analog  output  of  sufficient  magnitude  to  drive  an  x-y  recorder. 


Probe  Drive  Mechanism 


(ssnoS) 
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Figure  A.  Residual  leakage  field  profiles  for  rectangular  slots  of 

constant  width  (0.125")  and  varying  depth  after  dc  excitation 

a)  Typical  profiles 

b)  Peak-to-peak  magnitude  as  a function  of  depth 
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Typical  results  are  shown  in  Figures  3,  4,  5,  and  6.  Both  active 
and  residual  leakage  field  profiles  for  rectangular  slots  of  constant 
width  and  varying  depth  are  shown  in  Figures  3 and  4 respectively  to- 
gether with  a typical  lift-off  characteristic  and  plots  of  the  signals' 
peak-to-peak  magnitudes  versus  depth.  Corresponding  results  are  shown 
in  Figures  5 and  6 for  rectangular  slots  of  constant  depth  and  varying 
width.  Probe  "lift-off"  in  each  case  was  0.015". 

Note  particularly  the  differences  in  peak-to-peak  magnitude  between 
the  active  and  residual  leakage  fields,  as  well  as  the  slight  changes 
in  profile  shape. 

Care  was  taken  in  setting  up  the  experiment  to  ensure  that  a uniform 
magnetizing  current  distribution  was  obtained  over  the  cross  section  of 
the  specimen,  and  that  the  size  of  the  Hall  probe  was  sufficiently  small 
relative  to  the  size  of  the  specimen  that  a virtual  point-to-point  leak- 
age field  profile  was  obtained  by  scanning  over  the  specimen  surface. 


Dipole  Modeling 

Perhaps  the  most  widely  used  model  for  electromagnetic  field/defect 

interactions  is  that  developed  by  Zatsepin  and  Shcherbinln^^.  Initially 

proposed  for  the  prediction  of  residual  leakage  fields,  the  technique 

has  now  been  improved  and  extended  to  Include  both  dc  and  ac  excitation 
18—23 

conditions  . Leakage  fields  around  rectangular  slots  are  modeled  by 
assuming  a "uniform  magnetic  charge  distribution"  over  the  sides  of  the 
slot,  neglecting  material  nonlinearities  and  assuming  that  the  "charge" 
on  the  base  of  the  slot  does  not  contribute  to  the  leakage  field.  Under 
these  conditions  expressions  can  be  developed  for  the  x and  y components 
of  the  resulting  leakage  field  as  follows: 


Hx  * 2CTg[arctg 


h(Kfb) 


(x+b)  +y(y+h) 


- arctg 


h(x-b) 


-] 


(x-b)  4y(y+h) 


o In  + (y+h)^][(x-b)^  + y^] 

_ Xll  A ^ r\  fy 
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[(x+b)^  + y^][(x-b)2  + (y+h)^] 


(1) 


(2) 


where  the  slot  width  Is  2b  units  and  the  slot  depth  is  h units,  o^  is 

the  assumed  "surface  charge  density." 

Figures  7 and  8 show  plots  of  and  based  on  equations  1 and  2 

for  the  rectangular  slots  whose  actual  residual  leakage  field  profiles 

are  shown  in  Figures  4 and  6.  The  results  have  been  normalized  and  a 

12 

value  of  0 =1  assumed  in  each  case.  In  practice  , a values  are 

s s 

obtained  from  the  experimental  results  in  order  to  obtain  agreement 
between  the  theoretically  predicted  profiles  and  those  actually  measured. 
This  has  been  done  in  Figure  9 which  shows  a comparison  of  dipole  model 
and  experimental  results  for  a particular  slot  and  for  both  residual  and 
active  magnetization  conditions.  In  both  cases  the  experimental  peak 
values  have  been  substituted  into  equation  2 for  in  order  to  obtain 
the  plots. 

The  active  magnetization  case  has  also  been  Included  here  because 
the  dipole  model,  although  developed  initially  for  residual  magnetization 
conditions'^,  has  also  been  used  to  model  field/defect  interactions  with 
direct  current  excitation. 

Computer  code  generated  to  produce  the  dipole  model  leakage  field 
profiles  from  equations  1 and  2 is  given  as  Appendix  A.  The  program  is 
self-explanatory  with  regard  to  the  required  input  parameters  and  can 
be  used  to  predict  the  leakage  fields  around  rectangular  slots.  MAPA  is 
a plotting  routine  available  on  Colorado  State  University's  CDC  6400 


i 

I 
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I 
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computer  system. 


a)  Tangential  field  component,  H 

b)  Vertical  field  component,  H 

c)  Peak-to-peak  H magnitude  versus  depth 


Sli^n,.!  Ar,i.mude  (r.us.)  •npl^tude  ,gau»,i) 
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Finite  Element  Modeling 

Although  considerable  work  has  been  done  to  develop  closed  form 

mathematical  solutions  for  the  prediction  of  leakage  fields  around  defects 
2 A 23  26 

In  the  case  of  dc  and  ac  ’ excitation  conditions,  the  results  do  not 
lend  themselves  well  to  the  general  approach  to  defect  characterization 
called  for  in  the  introduction.  The  primary  reason  for  this  is  that, 
although  the  equations  governing  the  electromagnetic  field/defect  inter- 
actions are  readily  obtainable  in  a poissonian  partial  differential  equa- 
tion form,  the  nonlinear  magnetization  characteristics  of  ferromagnetic 
materials  and  complex  defect  geometries  found  in  practice  prevent  closed 
form  solutions  from  being  obtained  for  anything  but  the  simplest  condi- 
tions. 

Developments  in  numerical  analysis  technqiues  have  had  a significant 

impact  on  the  study  of  magnetic  flux  distributions  in  electrical  machine- 
27-31 

ry  and  the  principal  investigator  strongly  suggests  that  such  analysis 

techniques  could  have  a similar  impact  on  the  problem  of  defect  character- 
ization schemes  for  electromagnetic  methods  of  nondestructive  testing.  In 
the  case  of  a ferromagnetic  bar  carrying  a dc  magnetization  current  the 
vector  potential,  A in  the  region  is  related  to  the  current  density,  J 
and  reluctivity  (H/B) , v by  the  equation 


3 . 3A.  ^ 3 . 3A. 


(3) 


assuming  that  the  bar  is  at  the  center  of  an  unconventional  coordinate 
axis  system  and  carrying  current  in  the  positive  z direction.  The  finite 
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Figure  10.  Mesh  generation  for  the  finite  element  modeling  of  leakage 
fields  around  rectangular  slots  in  a ferromagnetic  bar. 
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2.  Minimizing  an  energy  functional.  In  this  case 

B 

W = //  [ /vbdb]dxdy  - //  JAdxdy  (4) 

R o R 

over  the  discretized  region  to  yield  a set  of  equations  for 

the  vector  potential  at  each  node  of  the  form 

A j = f(J,v,A  ,,  ^ , ,x,y)  (5) 

node  adjacent  nodes  ■' 

3.  Solving  the  nodal  vector  potential  equations  by  matrix 

Inversion  or  Iterative  techniques. 

The  finite  element  code  given  as  Appendix  B Is  based  directly  on 
32 

Anderson's  method  . Discussions  of  energy  functionals,  convergence, 
mesh  subdivision  and  B/H  curve  modeling  can  be  found  In  references  33 
to  37.  The  code  Is  self  explanatory  with  regard  to  the  required  Input 
data  and  consists  of  three  distinct  parts; 

1.  Mesh  generation  algorithm  (B.l) 

2.  Element  ordering  algorithm  (B.2) 

3.  Vector  potential  and  flux  density  calculation 
algorithm  (B.3) 

This  format  allows  maximum  flexibility  to  examine  a wide  variety  of 
geometries  although  the  three  parts  could  be  run  together  if  needed. 

In  examining  the  leakage  field  around  relatively  small  slots  in 
the  surface  of  a square  ferromagnetic  bar  it  is  possible  to  "zoom  In" 
on  the  region  of  interest  by  making  two  runs  with  the  computer  code. 
First,  a mesh  is  generated  as  in  Figure  10a)  for  the  square  bar  and  a 
surrounding  air  region  of  sufficient  size  that  the  vector  potential, 

A,  Is  zero  on  the  outer  boundary.  The  boundary  can  be  found  experi- 
mentally, using  a gaussmeter,  or  by  a trial-and-error  approach,  applying 
the  computer  code  several  times  until  no  further  change  In  A is  observed 


k 


8 

■ with  changes  in  outer  boundary  position.  After  calculating  the  vector 

potential  at  every  node  within  this  mesh,  the  "zoom  in"  region  around 
the  defect  is  chosen,  as  shown  in  Figure  10b) , with  the  vector  potentials 
on  the  region  perimeter  serving  as  boundary  conditions  for  the  second  run. 
Some  form  of  linear  interpolation  must  be  used  to  find  all  nodal  A values 
i for  this  second  run,  or  a hand-generated  mesh  could  be  used  similar  to 

I 

the  one  shown  in  Figure  11.  This  technique  was  used  by  Hwang  and  Lord 
for  a variety  of  defect  shapes  and,  based  on  the  results,  a defect 
characterization  scheme  was  proposed  for  electromagnetic  methods  of 
nondestructive  testing  relying  on  active  dc  magnetization  conditions'^ 

The  results  shown  in  Figure  13  were  obtained  using  the  mesh  of  Figure  12 
and  the  computer  program  listed  as  Appendix  B.3. 

Discussion 

In  order  to  develop  defect  characterization  schemes  for  electro- 
magnetic methods  of  nondestructive  testing,  a flexible  mathematical 
model  is  needed  to  predict  magnetic  field/defect  interactions  under  the 
constraints  of  nonlinear  operating  conditions  and  awkward  boundary 
shapes.  Although  the  dipole  model  predicts  the  shape  of  leakage  field 
profiles  around  rectangular  slots  in  ferromagnetic  materials,  the  under- 
lying linearizing  assumptions  associated  with  the  technique  prevent 
; its  use  as  the  basis  for  defect  characterization  schemes. 

Finite  element  analysis  techniques  have  been  applied  with  success 


J 


to  the  prediction  of  magnetic  fields  in  electrical  machines  and  this 


(normal)  component. 


technical  report  shows  how  they  can  be  used  to  predict  leakage  fields 
around  slots  in  ferromagnetic  materials  excited  with  a dc  magnetization 
current.  Certain  assumptions  do  have  to  be  made  with  regard  to  the 
location  of  the  zero  vector  potential  boundary,  discretization  and  B/H 
curve  model,  but  on  the  whole,  finite  element  techniques  are  sufficiently 
flexible  and  accurate  to  be  applied  to  the  defect  characterization  or 
"inverse"  problem. 

With  regard  to  extending  the  techniques  to  other  forms  of  excitation, 
the  following  comments  are  offered: 

Residual  leakage  fields  exist  around  defects  in  ferromagnetic 


materials  after  initial  magnetization  with  a dc  current,  causing 

regions  around  the  defect  to  behave  as  permanent  magnets  working 

in  the  second  quadrant  of  the  material's  B-H  loop  represented 

by  region  BC  in  Figure  1.  Some  work  has  been  done  in  studying 

flux  distributions  in  permanent  magnet  machines  using  numerical 
41-44 

analysis  techniques  which  could  be  applied  to  the  modeling 

of  leakage  field  profiles.  Experimental  studies  of  the  residual 
fields  within  slots  in  ferromagnetic  materials  are  being  carried 
out  by  the  authors  at  Colorado  State  University  to  define  in  detail 
the  permanent  magnet  behavior  of  the  material  in  the  vicinity  of 
the  slot.  The  results  will  provide  needed  modeling  data  for  the 
finite  element  program  which  In  this  case  provides  an  approximate 
solution  for  the  equation 
3 , 3A.  3 , 3A. 

is:  ^ ‘V  ■ \ 


where  J is  similar  to  the  current  density  term  of  equation  3, 

ID 

but  in  this  case  models  the  behavior  of  the  permanent  magnet 


material 


Eddy  currents  and  their  associated  fields  set  up  by  an  ac 
excitation  current,  react  with  a defect  to  produce  an  effec- 
tive change  in  the  impedance  of  a pick-up  coil  scanning  above 
the  surface  of  the  material  under  test.  In  this  case  the 
material  in  the  vicinity  of  the  defect  is  cycled  through  a 
complete  B-H  loop  (ABCDEA  in  Figure  1)  and  the  governing  equa- 
tion to  be  solved  by  numerical  techniques  is  of  the  form 


3x  dx  dy  3y  s 

where  is  the  imposed  current  density,  w is  the  angular 
frequency  and  a the  conductivity  of  the  material.  This  type 
of  equation  has  recently  been  solved  for  electrical  machine 
applications^^  using  finite  difference  and  finite  element 
analysis  and  such  techniques  could  also  be  used  to  study 
basic  eddy  current  phenomena  associated  with  electromagnetic 
methods  of  nondestructive  testing. 

Although  equations  3,  6 and  7 are  given  in  two  dimensional  form, 
51  52 

work  ’ would  indicate  that  three  dimensional  solutions  may  be 


(7) 


recent 

avail- 


able in  the  future. 
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Dipole  modeling  code  for  rectangular  slots. 
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*********‘f-  ********** 


PROGRAM  0IP0L3  CALCULATES  THE  MAGNETIC  FIELC  AROUND  A RECTANGULAR 
OEFECTt  USING  A STRIP  OIPOLt  MODEL. 

THE  REQUIRED  INPUT  CONSISTS  CF  THE  FCLLCWING  DATA. 

(SIGMA)  IS  A SCALING  CONSTANT  USED  TO  CESCRIBt  THE  LINEAR 

DENSITY  OF  THE  MAGNETIC  CHARGE  ALONG  THE  SURFACE  OF  THE 
DEFECT. 

(Y)  DESCRIBES  THE  LIFT-CFF  IN  INCREMENTS  CF  (INCR). 

(IDATA)  SPECIFItS  THt  TOTAL  NUMBER  CF  CEFtCTS  t>AMINtO. 

(MIN)  AND  (MAX)  ARE  THE  MINIMUM  AND  MAXIMUM  FCINTS  TO  BE 

SPECIFIED  TO  EITHER  SIDE  OF  THE  CENTER  CF  THE  DEFECT.  THE 
SPACING  BETWEEN  POINTS  IS  (INCRX). 

(HI)  AND  (W)  ARE  THE  CEPTH  AND  WIDTH  MEASURtMENTS  CF  THt 

DEFECT  (IN  UNITS  OF  (INCR)). 

A CROSS-SECTION  OF  THE  RECTANGULAR  SECT  IS  SHOWN  BELOW. 
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THE  OUTPUT  INCLUDES  THE  MAXIMUM  VALUES  FCUNC  FOR  THE  MAGNETIC 
FIELD  IN  BOTH  THE  X AND  Y DIRECTIONS  ((BX)  AND  ( BY ) , RESPtCT IV EL Y ) * 
THE  VALUES  USED  FOR  (SIGMA),  (MIN),  (MAX),  (Y),  AND  (INCR),  AND 
THE  PER  UNIT  VALUES  USED  FCR  EACH  DEFECT.  THIS  PROGRAM  ALSO 
PRODUCES  MICROFILM  PlOTS  CF  THE  MAGNETIC  FIELD  VALUES  (BX),  (BY), 
AND  OT)  AT  A SPECIFIED  lIFT-CFF.  (BT)  IN  THIS  CASE  IS  THt  TOTAc 
MAGNtTIC  FltLD,  OtFINtO  AS 

BT=(BX**2+BY**2) ♦♦O.S. 

PLOTS  OF  PEAK  TO  PEAK  (BY)  VALUES  VERSUS  DEFECT  DEPTH  ANC  WIDTH 
ARE  SHOWN  NEXT,  ACCOMPANIED  EY  TABlES  LISTING  THE  VALUES  USED  IN 
ALL  THt  PLOTS. 

«*««««»««*  *********** 


DIMENSION 

1 HX(5,241),HY(5,241),HXY(5,241),XX<2l1) 

2 ,H1(5) ,H(3) ,W(B) ,B(5) 

3 ,H2(241)  ,H3(241)  , H4  ( 2 '♦I ) , M5  (241  ) , H9  (24  1 ) 

4 ,HY4(5) ,HY5(5) ,XY1(3) , XY2 (5) 

* , MTIT(8) ,MTI2(8) ,MTI3 ( 1) 

REAL  INCR  , INCRX 

INTEGER  XT(8),YT(8),XXT(8),YTT(8),XTT(8) 

DATA 

1 XT/IOHX  OIST (IN) , 7*1H  /,YT/8*1H  / , X xT/ 1 0 HOEPTH 

2 ,YTT/10HP-P  ,7»1H  / ,YTT/10HWICTH  ,7*1H 


,7*1H 


uuu  oOu  uuu 


DATA  INPUT 


REAO(S,100)  SIGMA 

READ  (5,101)  Y,INCR,INCRX 

READ  (5,102)  IDATA,MIN,MAX 

READ  (5,103)  (Hl( I) ,M( I) ,1  = 1 .lOATA) 

——INITIALIZATION  OF  VARIABLES 

HY  3 = 0.0 
HX3=0.0 
MA=-MIN+MAX4-1 
N=0 

Y=Y*INCR 

DO  200  I=1,I0ATA 
HY4(I)=0.0 
hY5(I)=0.0 
H(I)  = 0.  0 
B(  i)=a.o 
200  CONTINUE 


calculatation  of  0x  and  3Y 
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c 


230 


2A0 

250 


DO  250  I=1,I0ATA 
H(I)=H1(I)*INCR 
3( I)= (W (I )/2)*lNCR 
DO  250  M=1,MA 

XX(M)= (FL0AT(M-1»MIN) ’INCKX) 

X1=XX(M) 

01= ( (X1+B(I) ••♦Z) ♦ (Y* ( Y>H(I) ) ) 

D2= ( (Xl-a(I) )**2) ♦(Y*( Y»H(  I) ) ) 

HX1=ATAN((H(I)*(X1+B(I)))/C1) 

HX2=ATAN( (H(I) * (X1-B(I ) ) )/02) 

HX( I,M)=2*SIGMA*(HX1-HX2) 

i 

IF(HX(  I,M) .GE.HX3.AND.HX (I,M).GT.0.0)  HX3  = HX(I,M)  \ 

i 

03=((X1>B(I))**2*(Y+H(I))**2)»((X1-B(I))**2+Y**2) 

D<*=  ( (Xl  + Bd)  ) •*2+Y**2)  • ( (X  1-a  (I ) )**2»(  Y*M  ( I) ) **2) 

HY (I,M)=SIGMA*AL0G(03/D4)  j 

j 

IF(HY(I,M) .GE.HY3.AN0.HY  (I,M).GT. 0.0)  HY3=HY(I,M)  i 

IF  (HY(I,M).LE.HY4(I).0R  . HY  (I  ,M) . LT. Q. 0 ) GO  TO  230  ’ 

HY^(I)=HY(I,M) 

XY1(I)=XI  1 

GO  TO  240  j 

CONTINUE  j 

IF(HY( I,M).Gt.HY5(I) .OR  . H Y ( I , F ) . G T . 0 . 0 ) GO  TO  240  ] 

HY5(I)=HY(I,M)  j 

XY2(I)=X1  I 

CONTINUE 
CONTINUE 


ooo  ooo  ooo 


c 

c-—— calculation  of  unit  values 
c 

00  2b0  I=1*I0ATA 
XNAX=A8S (HXi) 

YMAX=ABS (HY3) 

HY««(I)=HY<*(I>/YMAX 
HY5(I)=HY‘>(I)/YMAX 
DO  260  J = l,iiA 

HX(I, J)=HX(I,J) /XMAX 
HY(I, J)=HY(I,J» /YMAX 

HXY(I, J>=SQRT((HX(I, J»  >**2*(HY(I,J) »**2) 
260  CONTINUE 

PRINT  OUTPUT 


WRITE (6, 110) HX3,HY3 

WRITE  (6tll2)  SIGMAt MINtMAX 

WRITE  (6,113)  Y 

WRITE  (6,llii)  INCR 

DO  270  J=1,I0ATA 

WRITE  (6,111)  J,IOATrt 

WRITE  (6,115)  W(J) 

WRITE  (6,116)  HKJ) 

270  CONTINUE 


——PLOTS  AND  OATA 

267  DO  300  J=l,IOATA 
DO  300  M=l,f1A 

GO  TO  (277,278,279,200,  28A),J 

277  H2(M)=HX(J,M) 

GO  TO  300 

278  H3 (M)=HX ( J,N) 

GO  TO  300 

279  H<,(M)=HX(  J,M) 

GO  TO  300 

280  H5(M)=HX ( J,M) 

GO  TO  300 

284  H9(‘1)=HX(  J,M) 

300  CONTINUE 

READ  (5,109)  NTIT 

CALL  NAPA  (7,XX,H2,1,241«WL,NH, VL,VH,10HX  0IST(IN),1H 
CALL  NAPA  (7,XX,H3,1,241,WL,NH,VL,VH,10HX  0IST(IN),1H 
CALL  MAPA  (7,XX,H4,l,241,ML,MH,Vl,VH,10HX  CIST(IN),1) 
CALL  MAPA  (7,XX,H5,1,241,WL,WH,VL, VH,10HX  0IST(IN),1H 
call  MAPA  (7,XX,H9,1,241,WL,WH,VL,VH,10HX  0IST(IN),1H 


ii 

i 

1 


I 


i\ 


,HTIT,1) 

,NTIT,1) 

,MTIT,1) 

,MTIT,1) 

,MTiT,l) 


PRINT  FILM  PLOTS 


CALL  MAPM 
call  MAPM 
CALL  MAPM 
call  MAPM 
CALL  MAPM 
CALL  MAPM 
CALL  MAPM 
call  MAPM 
CALL  MAPM 
CALL  MAPM 
CALL  MAPM 
CALL  MAPM 


(6,XX,H2,l,2<4l,WL,WH,VL,VH,XT,YT,MTIT,l) 
(6,XX,H3,1,241,WL,WH,VL,VH,XT,YT,)<TIT,1) 
(6,XX,H4,1,2  41 , NL , WH,VL, VH,XT ,Y  T,MTIT,1) 
(o,XX,H5,1,2  41,WL,WH,  VL,  VH,XT,YT,r'TIT,  1) 
(b,XX,H9,l,241,WL,WH,VL,VH,XT,YT,MTIT,l) 
(1,XX,M9,1,241,WL,WH,VL,VH,XT,YT,PTIT,1) 
(2,XX,H2,1,2<41,Wl,WH,VL,VH,XT,YT,MTIT,1) 
(2,XX,H3,l,2  41,WL,WH,VL,VH,XT,YT,f'TIT,l) 
(2,XX,H4,1,2‘«1,WL,WH,VL,VH,XT,YT,MTIT,1) 
(2,XX,H5,1,241,WL,HH,VL,VH,XT,YT,MTIT,1) 
(2,XX,H9,1,241,Wl,WH,VL, VH,XT,YT,NTIT,1) 
(4,XX,H9,1,241,WL,WH,VL,VH,XT,YT,MTIT,1) 


J 


o o o 


C 


1 


N=N>1 

GO  TO  (350,351,362) ,N 

350  00  360  J=l,IOATA 
DO  360  M=1,MA 

360  HX ( J,M)=HY(J,M) 

GO  TO  267 

351  00  361  J=l,IOATA 
DO  361  M=1,MA 

361  HX  (J,ri)  = HXY{J,M) 

GO  TO  267 

352  CONTINUE 

00  400  I=ltIDATA 

HY4  (I  )=HY4*(I)  ♦ABS  (HY5(  I ) ) 

400  CONTINUE 

READ(5,109»  MT12 

CALL  MAPA ( 7,H1,HY4,1 ,5,WLt  Wh, VL, VH , 5HOEPT H , 3H P-P , MT I i , 1 ) 

CALL  NAPA (5,H1»HY4,1 ,5  » WL. WH, VL, VH , 5 hOE PT H , 3HP-P , PT I 2 ♦ 1 ) 

CALL  NAPN(5»H1,HY4,1,5 ,WL,WH,VL, VH,XXT, YTT,NTI2tl) 

450  CONTINUE 

READ  (5,109)  MTI3 

call  NAPA  (7,H,HY4,1,5,NL,WH,VL, VH,5HWI0TH,3HP-P,MTI3,1) 

CALL  NAPA  (5,W,HY4,1,5,WL,NH,VL, VH,5HWI0TH,3HP-P,MTI3,1) 

CALL  MAFM  ( 5 , W, H Y 4 , 1 , 5 , WL , WH , VL , VH , X T T , YT T , MT I 3, 1 ) 

500  STOP 

FORMAT  STATEMENTS 

100  FORMAT  (5X,F5.2» 

101  FORMAT  (5X,2F5,3,  F10,8» 

102  FORMAT  (5X,3I5> 

103  FORMAT  (5X,2F5.2) 

104  FORMAT  (110) 

109  FORMAT  (8A10) 

110  FORMAT  ( 1H1,5X,23HMAX  VALUES  FOUND  F CR  HX , F 1 0 . 4, 3X, 6 FFOR  HY,F10.4) 

111  FORMAT  (////, 5X,10HDATA  E N TR Y ,2 X , I 2 , 2 X , 2HOF  , I 5 ) 

112  FORMAT  (/,5X,7HSIGmA  = , F5 , 2 , 6 X, 1 5H X (MIN, /PAX.)  =,I5,1H/,I5( 

113  FORMAT  (/,5X,36H0ISTANCE  ABOVE  THE  SURFACE  CEFtCT  IS ,2X,F5, 3,2X,3H 
IIN.  ) 

114  FORMAT  (/,5X,33HTHE  DEFECT  DIMENSION,  IN  UMTS  OF  , 2 X , F5 . 3 , 2 X , 3HIN  . 
1) 

115  FORMAT  (/.‘♦2X,F5.2) 

116  FORMAT  (30X,11H— — — — — ,7X,11P— — — — ,3(/,40X,lH*,  7X,1H*)  ,2 

1X,F5,2,2(/,40X,1H*,7X,1H*),/,40X,  9 ) 

END 


Appendix  B 


Finite  element  modeling  code. 

1.  Mesh  generation 

2.  Element  ordering 

3.  Vector  potential  and  flux  density 
calculation 


o o o o o o o o o o o o o o o o n o n o o o o o o o o n r>  n n n o o n o 


& PROGRAM  MESMCINPUTfOUTPUTtPUNCH, FILMPLf TAP£5=INPU1. 

ITAPESzPUNCH) 

PROGRAM  MESH  GENERATES  A PLANAR  T WO- Cl MENSICNAl  TRI 
STRUCTURE  WITHIN  THE  SYMMETRICAL  BOUNDARY  SHAPE  SPE 
GIVEN  QOUNOARY  POINTS. 

IN  THE  ORDER  OF  READING,  THE  FCLLOWING  INPUT  DATA  I 
....INITNP  AND  INITEL  ARt  THE  INITIAL  NUMBERING  VAL 
FOR  THE  NODAL  POINTS  ANC  ELEMENTS  RESPECTIVELY, 
usually  1,1. 

....NOIVX  AND  NDIVY  ARE  THE  NUMBER  CF  ELEMENTAL 

DIVISIONS  ALONG  THE  X ANC  Y AXES  RESPECTIVELY. 
....XCOR  AND  YCOR  ARE  THE  (X,Y>  COuRCINATES  OF  THE 
EIGHT  DEFINING  BOUNDARY  POINTS  FCR  THE  GIVEN  SH 
....XMIN,  XMAX,  YMIN,  YMAX  ARE  THE  MINIMUM  AND  MAXI 
values  of  X AND  Y ASSOCIATED  WITH  THE  GIVEN  SHA 
(XMAX  - XMINI  SHOULG  BE  EQUAL  TO  (YMAX  - YMIN). 

OUTPUT  DATA  INCLUDES  THE  TOTAL  NUMBER  CF  NCCE  PCINT 
NUMBER  OF  ELEMENTS,  A LISTING  OF  THE  (X,Y)  COOROINA 
NODE  POINT  AND  THE  NODE  POINTS  DEFINING  EACH  ELEMEN 
WITH  A PLOT  OF  THE  MESH  STRUCTURE. 


DIMENSION  XCOR ( 8) , YCOR ( 3 ) , RN  ( 8 ) , XORD (1250  ) , YORC(125 

READ  AND  INITIALIZATION  OF  OATm 

REA0(5,1)  INITNP, INITEL 
WRI TE (6, 2) 

WRITE (o, 1) INI TNP, INITEL 

REA0(5, 3 )NOIVX, NDIVY 
WRITE(b,4) 

WRITE(6, 3)NOIVX, NDIVY 

REA  0(5, 5)  (XCOR(I) ,YCOR(I) ,1=1,8) 

WRITE (0,6) 

WRITE (6, 5)  (XCOR (I) , YCOR (I)  ,1  = 1,8) 

REA  0(5, 7 )XMIN, XMAX, YMIN, YMAX 
WRITE (6, e) 

WRITE(b,7)  XMIN, XMAX, YMIN, YMAX 

NUMNP=(N0IVX»1)»(N0I VY,1) 

NUMEL=2* (NOIVX) *( NDIVY ) 

WRITE (6, 13) NUMNP,NUMtL 

CALCULATE  NODAL  POINT  CCCRCINATES 

DX=NDIVX 
OX=l./OX 
DY=NQIVY 
0Y=1,/0Y 
C 

1EN0=N0IVV*>1 

JENO=NOIVX»i 

J1=0 


TAPE6=0UTPU7, 


***<nnt**** 

angular  element 

CIFIEO  BY  EIGHT 


S REQUIRED. 
UES 


AFE. 

MLM 

PE. 


S,  THE  TOTAL 
TES  FOR  EACH 

T,  TOGETHER 


0) ,NP(23eO,3) 


u u u 


C 


DO  250  I=l,IENO 

RY=R*OY 

C 

00  240  J=l,JtND 
r-( j-1) 
kX=R*DX 

RN{l)=1.0*(l.-ftX)  ♦d.-RY)*!  l.-2.*RX-2.*RY) 
RN(2)=4. 0*(RX)* (l.-RX»  *(l.-fiY) 

RN  (3)=-1.0*(RX)*(  l.-RY)  ♦ (1.-  i.  ♦R  ><•  2 . ♦R  Y ) 
RN(4)=4.0*(RX)*(RY)*(1.-RY) 
RN(5)=-l.»(RX)*(RY)*(3.-2. •RX-a.’RYJ 
RN(6)=4. *(RX) *(1.-RX) * (RY) 

RN( ?)=-!.♦ (l.-RX) ♦(RY) ♦ {1.>2.*RX-2.*RY) 

RN( 81=4. 0* (l.-RX) *(RY )♦ (l.-RY) 

C 

XORC( Jl)=0. 

YORD( J1)=0. 

C 

DO  230  K=l,8 

XORO<  Jl)  =XORO  ( Jl)  t-RNCK)  *XCCR  (K) 

YORO( J1)=Y0RD(J1) ♦RNCK) *yCCR (K) 

230  CONTINUE 
240  CONTINUE 
250  CONTINUE 

CALCULATE  NP  ARRAY 

12=0 
C 

DO  340  1=1, NOIVY 
C 

J1=(I-1)*(N0IVX4-1)4-1 

C 

00  330  J=l,NOIVX 

I2=I2>2 

11=12-1 

J2=  J14-1 

J4=Jlr(NOIVXTl) 

J3= J4*l 
C 

01=  (XORD (J3)-XORO (Jl))**2+  (YORC(J3)-YORD(J1))**2 
02= (XORO(J4)-XORD (J2) ) **2+ ( Y CRD ( J4) - YCRD ( J2 ) ) ♦*2 
C 

IF(C1.&T.02)  GO  TO  320 
C 

NPdl,  1)=J1 
NP(I1,2)=J3 
NP{ I1,3)=J4 
C 

NP(I2,1)=J1 
NP(I2,2)=J2 
NP(I2,i) =J3 
C 

GO  TO  329 

C 


320  CONTINUE 


c 

NPJ Iltl) =J1 
NP(I1,2)=J2 
KPdl,  3)=J4 
C 

NP( I2,l»  =J2 
NP(I2,2)=J3 
KP(I2,  3)=J4 
C 

329  CONTINUfc 
Jl=  J2 

330  CONTINUE 
340  CONTINUE 

00  1001  I=l,NUhtL 
00  1001  J=lt3 
NP(I,J)=NP(I, J)*INITNP-1 
1001  CONTINUE 
C 

C NRITE  STATEMENT 

C 

WRITE (6,9) 

IFAC=INITNP-1 
DO  420  I=1,NUMNP 
I1=I+IFAC 
C 

WRITE(b,10)  Il.XOROd)  , YCRC(I) 

WRITE(8,1h)  I1,X0RD(I)  .YOROd) 

420  CONTINUE 
C 

WRITE(6,11) 

IFAC=INITEL-1 
00  440  I=1,NUMEL 
I1  = H-IFAC 

WRITE(b, 12)11, (NP(I,J) , J=l,3) 
HRITE(8,15)Il,(NPd,J)  ,J=1,3) 

440  CONTINUE 
C 

C PLOT  ROUTINE 

C 

CALL  StT(0.0,1.0,0.0,1.0,X)'IN,XMAX,Y)'IN,YMAX,l) 
C 

DO  480  I=l,NUiEL 
C 

I1=NP(I ,1) 

I2=NP(I, 2) 

I3=NP(I, 3) 

C 

CALL  POINT(XORO(I1),VORC(II) ) 

CALL  VECT0R{X0R0d2)  ,YO«Od2)  ) 

CALL  VECTOR(XOROd3)  , Y0R0d3)  ) 

CALL  VtCTORlXORDdl)  ,YOROdl)  ) 

C 

480  CONTINUE 


u o u 


CALL  FRAME 
STOP 

FORMAT  STATEMENTS 
1 FORMAT(2I10) 


2 

FORMAT(20H1 

INITNP 

INITEL) 

3 

FORMAT<2I10) 

4 

FORMAT(20HO 

NOIVX 

NCI VY> 

5 

FORMAT(2F10.4l 

6 

FORMAT(20H0 

XCOR 

YCOR) 

7 

FORMAT(4F10.4) 

9 

FORMAT ( 40H0 

XMIN 

XMAX 

YMIN 

YMAX) 

9 

FORMAT( 30H0 

NP 

XCRO 

YORO) 

10 

F0RMAT( I10,2F10 .4) 

11 

FORMAT (40H0 

tLEM 

NP(1» 

N F ( 2 ) 

NP( 3) ) 

12  f ORMAT(AI10> 

13  FORMATC*  NUMNP=*,I5,*  NUMtL=*tI5) 

14  FORMAT(I10,20X,2F10.4) 

15  FORMAT( I10»10X,3I10) 

END 


ooooor>ooooooooof»oo«c»o 


PROGRAH  ELANP  LISTS  THE  ELEMENTS  AROLNO  EACH  NODE  POINT  INSIDE  THE 
GIVEN  BOUNDARY  SHAPE. 

INPUT  DATA  REQUIRED  FOR  THE  PRCGRAM  INCLUDES 

....NUMkL.NST.NEO  ARE  THE  TOTAL  NUMBER  OF  ELEMENTS 

GENERATED  FOR  THE  GIVEN  SHAPE  BY  THE  MESH  FROGRAN, 

THE  NUMBER  OF  THE  FIRST  NODE  POINT  INSIDE  THE 
GIVEN  BOUNDARY  AND  THE  NUMBER  OF  THE  LAST  NODE 
POINT  INSIDE  The  boundary  RESPECTIVELY. 

....A  total  LISTING  OF  THE  NODE  POINTS  CEFIMNG  EACH 

ELEMENT  AS  GENERATED  BY  LINES  420  TC  440  OF  FRCGRAM 
MESH. 

OUTPUT  DATA  CONSISTS  OF  A LISTING  OF  EACH  NOCE  POINT  BcTNctN  NST 
and  NEO  with  a count  OF  THE  NUMBER  OF  ELEMENTS  ARCUNC  EACH  NOCE 
POINT  ANC  THEIR  IDENTIFICATION. 

Those  node  points  between  nst  anc  nec  lying  cn  the  bcunoary  are 

lOtNTIFIABLE  BY  A ROW  OF  ZcRCS  AFTER  THt  NCCE  NUMBER. 


DIMENSION  NOP(20)«NP(2360t3l 
6 FORHAT(Ii0«10X,3I10> 

8 FCRMAKJIIO) 

REA0(5.8)  NUMELtNST.NEO 

REA 0(5, 6)  (N,NF(N,1),NP(N,2),NP(N,3) ,NN=1, NUMEL) 
M = 0 

DO  42  K=NST,NED 
ICONT=0 
DO  82  10=1,14 
NOP(IQ)=0 
82  CONTINUE 
KL=NST-2 
K<=K/KL 
K1=KK*KL«-1 
K2=KK*<L 

IF(K  .EG.  Kl)  GO  TO  43 
IF(K  ,EQ.  K2)  GO  TO  43 


DO 

50  KX=1, 

»2 

GO 

TO  (90,92),KX 

90 

MS  = 

K-KltM4 1 

IKK-1)*2» (NST-3) 

ME  = 

MS+3 

GO 

TO  94 

92 

MS= 

MS+2* (NST-3) 

ME  = 

MS*3 

94 

CONTINUE 

DO  44  I=MS, 

ME 

IF(NP(I,1) 

.EQ. 

K) 

GO 

TO 

78 

IF( NP(I ,2) 

.EQ. 

K) 

GO 

TO 

78 

IF(NP(I, 3) 

.EQ. 

K) 

GO 

TO 

78 

GO  TO  44 

78  iCONT  = ICONT«-l 
NOP (ICONT)=I 
44  CONTINUE 
50  CONTINUE 
M = M^1 

NSSSS=NST-5 
IF(M  .GT.  NSSSS)  M=0 
43  WRITE(6,eO)  K,ICONT, (NOP(N) ,N=1, 14) 
MRITE(d, 80)  K,1C0NT, (NQP(N) •N=l,14) 
42  CONTINUE 
so  F0RNAT(16I5) 

STOP 


ooo  o ooo  o o o o o o o o o o o o o o o o o o o o o o o o o o o o 


FROGRAH  POTEN* INPUT, OUTPUT, FILHPL#TAPE5=INPUT*TAPE6=CUrPUT) 

PROGRAM  POTEN  IS  A FIMTt  ELtMENT  ALGORITHM  WHICH  CALCULATES  THE 
VECTOR  potential  AT  EVERY  NCCE  OF  THE  REGICh  DISCRETIZED  INTO 
TRIANGULAR  ELEMENTS  BY  PROGRAM  MESH. 

POTEN  IS  BASED  ON  ANOERSONS  EQUATIONS  (REF. 32)  AND  REQUIRES  INPUT 
DATA  FROM  BOTH  THE  MtSH  ANC  ELANP  PROGRAMS  AS  WELL  AS  A MODEL  OF 
THE  B/H  CHARACTERISTIC  CF  THE  FERROMAGNETIC  MATERIAL  AND  THE 
BARS  CURRtNT  DENSITY. 

IN  GENERAL  THE  PROGRAM  HAS  TO  BE  WRITTEN  FCR  THE  PARTICULAR 
GEOMtTRY  UNDER  CONSIOERATI CN  AND  THE  RECUIREC  OUTPUT. 

INPUT  DATA  REQUIRED  INCLUDES  — 

NUMBGN,NUMOI V DEFINE  THE  SCANNING  PATH  FOR  THE  FLUX  DENSITY 

PLOTS 

NUMNP,NUMEL  ARE  THE  TOTAL  NUMBER  OF  NCOE  POINTS  AND  ELEMENTS 

RESPECTIVELY 

.....NAROS  AND  NAROE  CORRESPCNC  TO  NST  AND  NED  IN  THE  ElANP 
PROGRAM 

THE  LAST  THREE  RFAD  STATEMENTS  REFER  TC  THE  DATA  GENERATED 

SY  MESH  AND  ELANP  PROGRAMS 

RVPER  IS  THE  INVERSE  OF  THE  PERMEABILITY  CF  FREE  SPACE 

CURR  IS  THE  BAR  CURRENT  DENSITY 


DIMENSION  NP(116Q,3),CJ(1160)  ,RMU(  llEO)  ,SL(1160,E),MAT(llt)0) 

1 , A(630) , XORD (630),YORC(630) ,AREA (1160),BX(1160),8Y(lie0) 

2 ,B(llfaO)  ,IEDG(630  ),lETG(b30,8),BXC(50).BYC(50),3C(i>0) 

3 ,XX(50) 

INPUT  DATA 

READ  15,1)  NUMBGN,NUMr)lV 

READ (5, 2)  NUMNP,NUMEL, NAROS, NAROt 

READ(5,4)  (N,  A(N) ,X0R0(N) , YORC(N) ,NN=1,NUMNP) 

REA0(5,6)  (N,NF(N,1),NP(N,2),NP(N,3) ,NN=1,NUMEL) 

READ (5,8) (N,IEDG(N),(IET&(N, I ) , 1 = 1 . 8 ) , NN= NAROS,NARGt ) 

PVPER=7957Z4. 7155 
CURR=2.6016E>05 

SCALING  and  ELtMtNT  CALCULATIONS 

00  20  I=1,NUMNP 
XORO(I) =XORD(I) *. 0254 
YORD(I) =Y0R0(I)*. 0254 


20  CONTINUE 


o o o o 


i I 


f 

f 

i 

h 

! 

t 
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DO  38  I=1,NUMEL 
IP=NP(I,t) 

JP=NP(I,2) 

KP=NP(I ,3) 

XJ=XORO (JP)-XORD( IP) 
YJ=YORO(JP)-VORO( IP) 
XK=XORO{KP)-XORO(IP) 

Y<=YCRO( KP)-YORO(IP) 

ARE A(I)= (XJ*YK-XK*YJ) /2. 0 

Al=-( YK-YJ) /2,0 

Bl=YK/2. 0 

Cl=-YJ/2.0 

D1=(XK-X J)/2.0 

tl=-XK/2. 0 

Fl=XJ/2,  0 

SL (!♦!)= (A1*A1*01 ♦□!)/ area ( I ) 
SL(I.2)=(Bl»Bl*El*ei)/ARtA  <i ) 
SL(I.3)= (C1*C1+F1*F1)/AREA(I) 
SL(  I,4)  = (A1*BH-D1*E1)/AREA(  I) 
SL(I,5)=(Al*Cl*01»Fl)/ARtA(I) 
SL (I,6)=(BI*CI*E1*F1) /AREA( I) 
38  CONTINUE 


1 


•f 

i 

i 


INITIALIZATION  OF  AIR  ANO  IRON  REGIONS 
— CALCULATION  OF  BAR  CURRENT 


60 

DO  60  I=1»NUMEL 

RNU (I)=kWPER 
CJ(I)=0.0 

MAT ( I)=0 

CONTINUE 

C 

66 

IE3L=301 

CONTINUE 

661 

IE9R=IEBL^23 

DO  67  I=ItBLtIE3R 
CJ( I) =AREA(I) ♦CURR 

662 

e(i)si.65 

663 

RHU(I)32Q19.66 

666 

MAT II)*I 

67 

CONTINUE 

C 

c-.— 

IE8L  = IE9L*‘*8 

IF  (IEBL.LE.852) 

-CRACK  PORTION 

C 

68 

DO  68  1=791,796 
kMU( I)=RVPER 
CJ(I)=0. 0 

MAT (I)=0 

CONTINUE 

69 

00  69  1=839,862 
RMU<I)=RVPER 

CJ( i)=a. 0 

MAT (I)=0 

CONTINUE 

o o o o o r» 


SELECTION  OF  ACCELERATION  FACTOR 


ITER=0 
FAC=2.1 
ftFAR=0. 07 

70  CONTINUE 
IFLA6=0 
ACK=0.0 
ITER=ITER*1 

IF  <ITER.EQ.6J  GO  TO  71 

FAC=FAC-Q.2 

GO  TO  72 

71  FAC=t.O 

72  CONTINUE 

eEGINING  OF  ITERATION 

00  326  KOA=1.40 
DO  42  K=NAROS.NAROE 
IFdEOGCK)  >201,42,201 
201  T1=T2=T3=0.0 
ME=ItO&>K) 

DO  44  JJ=1,ME 
I=IETG(K, JJ) 

IF  (NP(I,1)-K)  77,  78,77 

77  IF(NP(I,2)-K)  79,80,79 

79  IF(NP(I, 3)-K>  44,82,44 

78  IP=1 
JP=4 
KP=5 

I1=NP(I,2) 

I2=NP(I,3) 

GO  TO  84 

80  IP=2 
JP  = b 
KP=4 

I1=NP(I,3) 
i2=NP(I ,1) 

GO  TO  84 
82  iP=3 
JP=5 
KP=6 

I1=NPCI,1) 

I2=NP(I,2) 

84  T1=T1+CJ(I) 

12=T2^RHU(I)*SL(Z«1P) 

T33T34>RHU(I)*(Aai>*SL(I,  JP)  ♦A(I2)  *SL(I  ,KP)  I 
44  CONTINUE 
Tl*Tl/3.0 
AZ*(T1-T3»/T2 
ANEH=A(K)>FAC*«AZ-A(K) ) 

IF  <A3S(A(<>-ANEH) .Lt.  0.000  1 ) GO  TO  41 
IFLAG=1 

ACK=ABS (A(K)-ANEH) 

41  CONTINUE 
AIK>=ANE  R 
42  CONTINUE 

KNAs(KOA/2)*2 

IFfKHA  .HE.  KOA)  60  TO  326 


OOOO  ft  ftftft  Oftft 


CALCULATION  OF  NEW  FLUX  DENSITY  IN  THE  MATERIAL 


00  52  I=l,NUMeL 
IF(MAT(I)-0)  203,52,203 

203  IP=NP(I,1» 

JP=NP(I,2) 

KP=NP(I,3) 

BX (I)=(A( IP» ♦(XORO(KP) -XORC( UP) ) ♦Al JF) *(XCftClIF)-XCRC (KP) )*A(KP)* ( 
IXORDC JP>-XORO(IP» ) )/(2.*AREA(I)  ) 

BY  (I)  = (AdPJ*  (YORD(KPJ-YORC(  JP)  ) ♦A  (JP)*  (YCRC  (IP)-YCRC(<Pn  ♦A  (KP»  ♦ ( 
2YORO(JP)-YORO(IP) ) ) / ( 2 . * ARE  A ( I » ) 

B(I»=SQRT(aX(I)*BX(I>>dY(I)»eY(I)  ) 

CALCULATION  OF  NEW  PERMEABILITY  USING  MCCElEC  6/H  CURVE 

IF(B(I»-1.3»  91,91,204 

204  IF(a(I)-1.9)  93,93,205 

205  IF(a(I»-2.2)  95,95,206 

206  IF(B(I)-2.5»  97,97,207 

207  RR1=RVPER 
GO  TO  99 

91  RR1  = 95.0«-240.5576923*B(I) 

GO  TO  99 

93  RRl  = a(I) *(89623.8  2-B(I»* (7  AO  98.4 562-20369, 373*d(I>) >-35628.  365 
GO  TO  99 

95  RRl=6375,0rbl628,78767»(e(I)-1.9) 

GO  TO  99 

9 7 RR1=2536  3.6363»17  945««.545  7 * (3  (I) -2. 2) 

99  RMU  (I)=RMU(I>-*-RFAR*<RRl-RMU(I>> 

52  CONTINUE 
326  CONTINUE 

IF  (ITER.EQ.6)  GO  TO  353 
IF  (IFLAG.EQ.l)  GO  TO  70 
350  CONTINUE 

WRITE  (6,101)  ITER,ACK 


calculation  of  new  flux  density  in  air  portion 

DO  3500  I=1,NUMEL 
IF  (MAT(I)-a)  3500,400,3500 
400  1P=NP(I,1) 

JP=NP(I,2) 

KP=NP(I,3) 

OX (I) =( A(IP>* (XORO(KP)-XORC ( JP) ) ♦A (JP) * ( X ORC ( IP ) -XCRC ( KP ) ) ♦A(KP) * ( 
1X0R0( JP)-X0R0(IP>  >1/(2. *ARE AID) 

BY(I)=( A(IP) •(YORO(KP)-YCRD( JP) ) ^A ( JF ) * ( YC RC ( IF )- YOR C ( KP) )^A(KP) • ( 
2YORO(JP)-YORD<IP) ) )/(2.*AREA  (I) ) 

3500  CONTINUE 


ooo  ooo  o n o non 


CALCULATIOM  OF  NAGNCTIC  FIULO  ABOVE  BAR 
ICONTsQ 

NUHEN0=NUMBGN*NUM0IV*2-1 
00  510  I=NUHBGNtNUHEN0«2 
IC0NT=ICCNTfl 

BXC(IC0NT)=(BX(I)*BX(I+l))/2. 

BYC (IC0NT)= (BY(I) ♦BY(I+l))/2, 

BC(IC0NT»=SQRT(BXC(IC0NT)*exC(IC0NT)  ♦BYC(ICCM)*BYC(  ICCNT) ) 
XX (ICONT)=XORD (ICONT) 

COMINUt 

OUTPUT  DATA  AND  PLOTS 


CALL 

MAPA 

(7, XX,aXC,l,24,WL,WH,VL,VH 

,1H 

• IH 

tlH 

.1) 

CALL 

MAPA 

(5,XX,BXC,1,24,WL, WH,VL, VH 

,1H 

, IH 

.IH 

.1) 

CALL 

MAPM 

(5,XX,BXC,1,24,WL,WH,VL, VH 

tlH 

, IH 

,1H 

.1) 

call 

MAPA 

(7,XX,BYC,1,24,WL,WH,VL,VH 

,1H 

,1H 

.IH 

.1) 

CALL 

MAPA 

(5,XX,BYC,1,24,WL,WH,VL, VH 

, IH 

.IH 

.1) 

CALL 

MAPM 

(5,XX,8YC,l,2  4,WL,WH,V  L, VH 

,1H 

,1H 

.IH 

.1) 

call 

MAPA 

(7,XX,BC  ,i,24,WL,WH,VL, VH 

,1H 

»1H 

.IH 

.1) 

CALL 

MAPA 

(5,XX,BC  ,1,24,ML,WH,VL,VH 

,1H 

, IH 

.IH 

.1) 

CALL 

MAPM 

(5,XX,BC  ,1,24,WL  , )>H,VL,  VH 

,1H 

f IH 

.IH 

,1) 

WRITE  STATEMENT 
WRITE  (6,10» 

WRITE  (6,12)  (HAT (I) ,I=l,NUMtL,2  ) 
WRITE  (6,10) 

NUMB=1 

00  550  J=l,25 
NUM£=NUM3t’12 

WRITE  (6,11)  (A(I),I=NUMB, hUME) 
NUMB=NUM6^25 
CONTINUE 
NUMB=14 
CO  560  J=l,25 
NUME=NUM3>11 

WRITE  (6,14)  (A(I ) ,I=NUMB,NUM£) 
NUMe=NUMd^25 
CONTINUE 
STOP 

FORMAT  STATEMENT 

I FORMAT  (2110) 

2  FORMAT(4I10) 

4 FORMAT(Il0,10X,3F10.4) 
b FORMATd  10,10X,3I10) 

8 FORMATdOIS) 

10  FORMAT  (IHl) 

II  FORMAT  (1X,13£10.3) 

12  FORMAT  (2414) 

14  FORMAT  (3X,12(1E11.4) ) 

16  FORMAT  (lX,I10,10X,3Ei6.6) 

101  FORMAT  (1X,I10,E15.6) 

102  F0RMAT(5X,b<lX,F8.5) ) 

106  rORMAT(///) 

END 


